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We present full accounts of a method to extract nucleon-nucleon (NN) potentials from 
the Bethe-Salpter amplitude in lattice QCD. The method is applied to two nucleons on 
the lattice with quenched QCD simulations. By disentangling the mixing between the S- 



state and the D-state, we obtain central and tensor potentials in the leading order of the 
velocity expansion of the non-local AW potential. The spatial structure and the quark mass 
dependence of the potentials are analyzed in detail. 
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§1. Introduction 
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The origin of the nuclear force is one of the major unsolved problems in particle 

and nuclear physics even after the establishment of the quantum chromodynamics 

(QCD). Although the nuclear force is still not well-understood theoretically, a large 

number of proton-proton and neutron-proton scattering data as well as deuteron 

lO ' properties have been accumulated and summarized e.g. in the Nijmegen database.^ 1 

00 

. pion production threshold together with the deuteron properties, the notion of the 

\ NN potential (either in the coordinate space or in the momentum space) turns 

■ out to be very useful: 2 ' it can be determined phenomenologically to reproduce the 

scattering phase shifts and bound state properties either through the Schrodinger 
equation for the NN wave function or through the Lippmann-Schwinger equation 
for the NN T-matrix. Once the potential is determined, it can be used to study 
systems with more than 2 nucleons by using various many-body techniques. 

Phenomenological NN potentials which can fit the NN data precisely (e.g. 
more than 2000 data points with x 2 /dof ~ 1) at Ti a b < 300 MeV are called the high- 
precision iViV potentials: They include the potentials such as CD-Bonn, 3> Argonne 
ui8,® and Nijm I, Nijm II and Reid93.^ Also systematic low energy construction of 
the nuclear force on the basis of the chiral perturbation theory is being developed.®'^ 
The phenomenological NN potentials in the coordinate space are known to 
reflect some characteristic features of the NN interaction at different length scalesP 
(i) The long range part of the nuclear force (the relative distance r > 2 fm) is 
dominated by the one-pion exchange introduced by YukawaP 1 Because of the 
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pion's Nambu-Goldstone character, it couples to the spin-isospin density of the 
nucleoli and hence leads to a strong spin-isospin dependent force, namely the 
tensor force. 

(ii) The medium range part (1 fm < r < 2 fm) receives significant contributions 
from the exchange of two-pions (tttt) and heavy mesons (p, u, and a). In 
particular, the spin-isospin independent attraction of about 50 - 100 MeV in 
this region plays an essential role for the binding of atomic nuclei. 

(iii) The short range part (r < 1 fm) is best described by a strong repulsive core as 
originally introduced by JastrowP Such a short range repulsion is important 
for the stability of atomic nuclei against collapse, for determining the maximum 
mass of neutron stars, and for igniting the Type II supernova explosions.^ 1 

(iv) There is also a strong attractive spin-orbit force in the isospin 1 channel at 
medium and short distances. This leads to the 3 P2 neutron pairing in neutron 
matter and hence the neutron superfluidity inside neutron starsP^ 1 

A repulsive core surrounded by an attractive well is in fact a common feature of 
the "effective" potential between composite particles. The Lenard-Jones potential 
between neutral atoms or molecules is a well-known example in atomic physics. The 
potential between 4 He nuclei is a typical example in nuclear physics. The origin of the 
repulsive cores in these examples is known to be the Pauli exclusion among electrons 
or among nucleons. The same idea, however, is not directly applicable to the NN 
potential, because the quark has not only spin and flavor but also color which allows 
six quarks to occupy the same state without violating the Pauli principle. To account 
for the repulsive core of the NN force, therefore, various ideas have been proposed as 
summarized in Ref. Hip : exchange of the neutral u mesonp^* exchange of non-linear 
pion field,^ and a combination of the Pauli principle with the one-gluon-exchange 
between quarksPS Despite all these efforts, convincing account of the nuclear force 
has not yet been obtained. 

In this situation, it is highly desirable to study the NN interactions from the 
first principle lattice QCD simulations. A theoretical framework suitable for such 
purpose was first proposed by LiischerP^ For two hadrons in a finite box with the 
size L x L x L in periodic boundary conditions, an exact relation between the energy 
spectra in the box and the elastic scattering phase shift at these energies was derived: 
If the range of the hadron interaction R is sufficiently smaller than the size of the box 
R < L/2, the behavior of the two-particle Bethe-Salpeter (BS) wave function tp(r) 
in the interval R < \r\ < L/2 under the periodic boundary conditions has sufficient 
information to relate the phase shift and the two-particle spectrum. 

Liischer's method bypasses the difficulty to treat the real-time scattering process 
on the Euclidean latticeo Furthermore, it utilizes the finiteness of the lattice box 
effectively to extract the information of the on-shell scattering matrix and the phase 
shift. This approach has been applied to extract the NN scattering lengths in the 

*' There are several studies of the NN interactions on the lattice, which do not rely on Liischer's 
method. One uses the Born-Oppenheimer picture, i.e., if one of the three quarks inside the baryon 
is infinitely heavy, one may define the potential between baryons a la Born-Oppenheimer.^ This 
is, however, not applicable to the nucleons with light quarks. The other employs the strong coupling 
limit, which has been proposed quite recently" 
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quenched QCD simulations® and in the (2+l)-flavor QCD simulations with the 
mixed action.^ 1 

Recently, the present authors proposed a closely related but an alternative ap- 
proach to the AW interactions from lattice QCDP^*^ The starting point is the 
same BS wave function il){r) as discussed in Ref. [T5j) . Instead of looking at the 
wave function outside the range of the interaction, we consider the internal region 
|r| < R and define the energy-independent non-local potential U(r, r') from ip(r) so 
that it obeys the Schrodinger type equation in a finite box. Since U(r, r') for strong 
interaction is localized in its spatial coordinates due to confinement of quarks and 
gluons, the potential receives finite volume effect only weakly in a large box. There- 
fore, once U is determined and is appropriately extrapolated to L — > oo, one may 
simply use the Schrodinger equation in the infinite space to calculate the scattering 
phase shifts and bound state spectra to compare with experimental data. Further 
advantage of utilizing the potential is that it would be a smooth function of the 
quark masses so that it is relatively easy to handle on the lattice. This is in sharp 
contrast to the scattering length which shows a singular behavior around the quark 
mass corresponding to the formation of the AW bound stateF^l 

Since we consider the non- asymptotic region (|r| < R) of the wave function, the 
resultant potential U and the T-matrix are off-shell. Therefore, they depend on the 
nucleon interpolating operator adopted to define the BS wave function. This is in a 
sense an advantage, since one can establish a one-to-one correspondence between the 
nucleon interpolating operator and the AW potential in QCD, which is not attainable 
in phenomenological NN potentials. It also implies that the AW potential on the 
lattice and the phenomenological NN potentials are equivalent only in the sense that 
they give the same observables, so that the comparison of their spatial structures 
should be made only qualitatively. 

The purpose of this paper is twofold: First, we will present a theoretical foun- 
dation of our method to extract the NN potentials from lattice QCD. Then, we 
will give a full account of the application of the method to the quenched lattice 
QCD simulations. Once our method in lattice QCD is proved to work in the AW 
system, it will have various applications not only to nuclear many-body problems 
but also to hyperon-nucleon, hyperon-hyperon and three-nucleon interactions which 
have much less experimental information than the AW systems. A first attempt to 
the hyperon-nucleon potential has been already reported in Ref. and more on 
hyperons will appear in the future publications. 

This paper is organized as follows. In £J21 we illustrate the derivation of the 
two-body and many-body potentials from the wave function in quantum mechanics. 
In 331 the idea in the previous section is generalized to the interaction of composite 
particles in field theory. In £JH we classify the general structure of the AW potential 
in the velocity expansion and show the procedure to determine each term. In ^3 the 
method to determine the AW potential from the lattice QCD data is discussed in 



*' Similar situation is well-studied in connection with the BEC-BCS crossover in cold fermionic 
atomsp^ where the external magnetic field plays a role of the quark mass in QCD. For seminal 
suggestion on the rapid quark-mass dependence of the NN scattering length, see Ref. I23[l. 
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detail for the effective central potential at low energy. We also discuss the method 
to extract the tensor potential in our approach. In £}6l NN potentials obtained 
from the quenched lattice QCD simulations are presented. Section [7] is devoted to 
summary and concluding remarks. In Appendix [Aj a field-theoretical derivation of 
the asymptotic BS wave function at large distance is presented. In Appendix |Bl 
the way to make general decomposition of the NN potential (the Okubo-Marshak 
decomposition® ) is reviewed. In Appendix [Cj matrix elements of the general NN 
potential are presented. In Appendix D, heat-kernel representation of the Green's 
function is presented. 



§2. Non-local potential in quantum mechanics 



2.1. Two-body force 

To show the basic concept of the non-local potential in a finite box with the 
size L x L x L, we start with a non-relativistic two-body problem described by the 
stationary Schrodinger equation: 

(V 2 + k 2 n )^ n {r) =2fxJ U(r, r>«(r')dV ', (2-1) 

where r is the relative coordinate of the two spinless and non-relativistic particles, 
and k n is related to the discrete energy eigenvalues E n = k\j (2/x) (n = 0, 1, 2, • • • ) 
with \i being the reduced mass. The wave function obeys the periodic boundary 
condition. The non-local potential U(r, r / j*l is assumed to be energy-independent 
and Hermitian, U*(r',r) = U(r,r'), so that the discrete energy eigenvalues E n are 
real and corresponding eigenfunctions can be made orthonormal. For the scattering 
states (bound states) in the infinite volume, we have E(L — > oo) > (E(L —> oo) < 
0). On the other hand, negative E n (L) in the finite volume does not necessarily 
imply the existence of the bound state at L — > oo. 

We consider the potential whose spatial extension is sufficiently small in the sense 
that U(r, r') is exponentially suppressed for \r'\} > R with R being smaller than 
L/2. We define the "inner region" by f2 in = {r G L 3 | \r\ < R}. Then, the wave 
function in the "outer region" I? ut = L ?J — Q\ n satisfies the Helmholtz equation, 
(V 2 + /c 2 )V ; n(^) = 0, with the periodic boundary condition. 

Let us consider the following inverse problem: Suppose we have no information 
about U except that it is smooth and short ranged, while we know linearly indepen- 
dent wave functions ipn(^) and associated energy E n = k\j(2n) in a finite box for 



n < nA**}\ Now, we introduce the following function: 



K n {r) = -^(V 2 + k 2 n )Mr) = (r\(E n - H )\n), (2-2) 
where Hq is the non-relativistic kinetic energy operator satisfying {r\Ho\n) = 



*' Here we use the standard term "non-local" in the sense that U(r,r') cannot be written as 
V(r)S(r-r'). 

**' This is more luxurious situation than the usual inverse scattering problem where only the 
scattering phase shifts in the outer region are available. 



Lattice Nuclear Force 



5 



2^V 2 'f/ ; n( r ')- Since (V 2 + /c 2 ) removes the non-interacting part of the wave func- 
tion, K n {r) is non- vanishing only in the inner region J7j n irrespective of the sign of 
h 2 

By taking into account the fact that i^ n { r ) = (r\n) may not be orthonormal, we 
introduce the norm kernel M nn > = (n\n') = J (f 3 7"0* (r)ip n i(r), so that the projection 
operator to the space spanned by the wave functions with n < n c reads P(n c ) = 
X^raV l n )^Cn'( n/ | = T^n Pn- Then, an energy-independent and non-local potential 
can be defined as 



U(r,r') = (r\ 



H())Pn 



J2 K n(r)Kn'r n >(r'), (2-3) 



which leads to the Schrodinger equation Eq. (|2-ip for ip n <n c (. r )- If we apply a unitary 
transformation A to the wave function, ip — >• ijj' = Aip, the non-local potential is 
modified as U — > V = AUA^. Such unitary transformation does not affect the 
observables, while it changes the spatial structure of the wave function and the non- 
local potential. 

If E n are all real and N nn > = 8 nn >, the potential U = Y^i(E n — Ho)P n becomes 
a hermitian operator (n\U\n')* = (n'\U\n) in the subspace n < n c . Otherwise, 
the hermiticity is not obvious and should be checked case by case. In field theory 
discussed later, ip n (r) corresponds to the equal-time Bethe-Salpeter amplitude in a 
finite box and E Uc corresponds to the threshold energy Eth of inelastic channels. 

In practice, the potential defined in Eq. (|2-3|) has limited use, because the number 
of states satisfying the condition E < £th is not generally large for lattice QCD in a 
finite box. This problem can be evaded when we focus on the low-energy scattering 
with E sufficiently smaller than the intrinsic scale of the system or the scale of the 
non-locality of the potential. In such a case, the velocity expansion of U(r,r') in 
terms of its non-locality is usefulP^ 1 For example, a spin-independent potential with 
hermiticity, rotational invariance, parity symmetry, and time-reversal invariance can 
be expanded as 



U(r,r') = V(r,v)S(r-r'), 

V(r,v) = V (r) + ^{V v2 (r),v 2 } + V e2 (r)L 2 + 



(24) 
(2-5) 



where v = p/fi and L = r x p with p = —iSJ. Each coefficient of the expansion is 
the local potential and can be determined successively by the wave functions at low 
energies: For example, if we have five wave functions corresponding to £Vi=o,i,2,3,4, 
we obtain 



(E n - H )i> n (r) 



V (r) + ^{VAr),v 2 } + V e2 (r)L 2 



(2-6) 



Pretending that V v2 {r) and (-§^,) n V v2 (r) are independent of each other, Eq. (12-6)) for 
n = 0, • • ■ ,4 can be solved algebraically to obtain Vo(r), V v i(r), -§^.V v2 {r), (-§^) 2 V v2 (r) 
and Vp (r) . Hermiticity of the potential can be checked by the consistency among 



6 



S. Aoki, T. Hatsuda and N. Ishii 



the local potentials thus determined. Stability of the potentials against the number 
of wave functions introduced can be also checked. 

An advantage of defining the potential from the wave functions in the "inner re- 
gion" is that the effect of the periodic boundary condition is exponentially suppressed 
for finite range interactions: Then one can first make appropriate extrapolation of 
U(r,r') or V(r,v) to L — > oo, and then solve the Schrodinger equation using the ex- 
trapolated potential to calculate the observables such as the phase shifts and binding 
energies in the infinite volume0 This is in contrast to the approach by LiischerSS 
in which the wave functions in the "outer region" suffering from the boundary con- 
ditions is ingeniously utilized to probe the scattering observables. Apparently, the 
two approaches are the opposite sides of a same coin. 

2.2. Many-body forces 

For the interactions among composite particles, there are in principle many- 
body forces which take place in the system composed of more than two particles. 
The well-known example in nuclear physics is the Fujita-Miyazawa type three-body 
force acting among three nucleons.^'GIil 1 It is phenomenologically important for the 
extra binding of light nucleP^ and for the extra repulsion in high density matter^ 
and in elastic nucleus-nucleus scatterings.^ 

The method to define the two-body potential from the relative wave function 
discussed above can be generalized to the many-body forces. Let us illustrate the 
procedure by considering the three-body system of spinless and distinguishable par- 
ticles with equal mass m. We consider the local potentials for both two-body and 
three-body forces just for simplicity. In the rest frame of the three-body system, we 
have 



(E n - H 0r - H 0p )tp n (r, p) 



^2v 2 (xi,Xj) + ^3(011,052,053) 

i>j 



Mr,p), (2-7) 



where r(= x\ — x 2 ) and p(= 033 — (x\ + x 2 )/2) are the Jocobi coordinates. Hq t = 
— V^/(2/x r ) and Hq p = — V 2 p /(2p, p ) are the kinetic energy operator with p r {= m/2) 
and p p {= 2m/3) being the reduced masses. E n is the total energy of the three-body 
system at rest. Because of the translational invariance, the two-body potential V 2 
and the three-body potential V3 are the functions of r and p. 

If we know the wave function and the total energy on the left-hand side of 
Eq. (|2-7[) . the three-body potential can be determined by the following procedure. 
We first consider the situation, \p\ 3> R S> \r\, where V^a^,^), V^i^i,^) and 
1/3(031,032,053) are vanishingly small because of the assumed short-range nature of 
the potentials. Then, V 2 {x\,X2) can be determined by changing r within the range 
R > \r\. One can carry out similar procedure to determine V 2 (x 2 , 033) and V 2 (xi, 033). 
Alternatively, one may determine V 2 from the genuine two-body system. 

Once all the two-body potentials are determined, V3 can be extracted from 
the wave function in the range, R > \r\ and R > \p\, through the three-body 

*' Strictly speaking, the local potentials with higher derivatives must be treated as perturbation 
to keep the Schrodinger equation as a second order differential equation. 
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equation Eq. ()2-7j) . It is important to note that the three-body potential is always 
obtained together with the two-body potential: they are closely tied through the 
wave function. If one makes the unitary transformation of the wave function, both 
Vz and V3 are changed simultaneously. 

The above procedure can be formally generalized to the non-local potentials and 
to the iV(> 3)-particle systems with different masses and internal degrees of freedom. 

§3. Non-local potential in field theory for spin 1/2 particles 

3.1. Bethe-Salpeter wave function 

In field theory, the best analogue of the two-particle wave function is the equal- 
time Bethe-Salpeter (BS) amplitude, so that we use the term "BS wave function" 
throughout this paper. Let us consider the following BS wave function for the 6- 
quark state with total energy W and the total three-momentum P = in a finite 
box L 3 , 

&ap(r,t) = (0\n p (y,t)p a (x,t)\B = 2;W,P = 0) =^ a p{r)e~ iW \ (3-1) 

where the relative coordinate is denoted as r = x — y. The local composite operators 
for the proton and the neutron are denoted by p a (x,t) and n@(y,t) with spinor 
indices a and (3. The QCD vacuum is denoted by |0), while the state \B = 2; W, P = 
0) is a QCD eigenstate with baryon number 2 and with the same quantum numbers 
as the pn system. One should keep in mind that \B = 2; W, P = 0) is not a simple 
superposition of a product state |p) <g> |n), since there are complicated exchanges of 
quarks and gluons between the two composite particles. The stationary BS wave 
function ip(r) may be regarded as a probability amplitude in \B = 2; W, P = 0) to 
have "neutron-like" three-quarks located at point y and "proton-like" three-quarks 
located at point x. 

The spatial extent of the NN interaction in QCD is short ranged and is expo- 
nentially suppressed beyond the distance R > 2 fm. Therefore, the spatial part of 
the BS wave function in the "outer region" (r > R) satisfies the Helmholtz equa- 
tion, ((W/2) 2 - V 2 + m 2 N )^ a p(r) = -(V 2 + k 2 )ip a f3(r) = 0, up to an exponentially 
small correction. Here the "asymptotic momentum" k is related to the total en- 
ergy W through the relation, W = 2^J k 2 + m 2 N . To make a formal resemblance 
with the non-relativistic case, we introduce the "effective center of mass energy", 
E = k 2 /(2fi) = k 2 /mM^^ As shown in Appendix A, using the unitarity of the 
5-matrix, we can show that the asymptotic behaviour of the BS wave function at 
large r is identical to that of the scattering wave in the quantum mechanics, with 
the identification that the phase of the S-matrix is the scattering phase shift of the 
BS wave function. 

Now, we apply the same logic as the quantum mechanical case in §2.11 The 
threshold of the pion production _E t h — m n is chosen to be E nc . Namely, (V 2 + 
k 2 )il> a ^E{ r ) is a function which has a support only in the inner region as long as E 
stays below the threshold. Thus we can define the short-ranged non-local potential 
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as 

(E - H )i/>afiM r ) = J u *P;js(r, r')^ 7 5,£(r')dV , (3-2) 
U a p- jS (r, r') = K a p, E (rWEE^5,E> (O (3-3) 

E,E' 

= V a p n s(r,v)S(r-r'), (3-4) 

where E = k 2 jm^ and Hq = — V 2 /mjv- By construction, the solution of Eq. (j3-2p 
with U a /3 n s(r,r') extrapolated to L — > oo reproduces the correct BS wave function 
in the asymptotic region, and hence the phase shifts and binding energies of the 
two-nucleon system. 

The Schrodinger type equation with the non-local potential similar to Eq. (|3-2p 
has been derived for bosons on the basis of a diagrammatic method in Refs. [T5]) and 
[32]) . A slight difference is that our non-local potential has no explicit -E-dependence 
by construction as seen in Eq. (|3-3p . 

3.2. Interpolating operators 

In Eq. (|3-ip . simplest interpolating operators for the neutron and the proton 
written in terms of the up-quark field u(x) and the down-quark fields d{x) would be 

np{y) = e abc {u a {y)C^d b {y)) d c p(y), (3-5) 

Pa(x) = £ abc (u a {x)C^d b (x)) U ca (x), (3-6) 

where x = (x,t), y = (y,t) and the color indices are denoted by a, b and c. The 
charge conjugation matrix in the spinor space is denoted by C. 

As shown in Appendix A, local operators such as given in Eqs. (|3-5p and (|3-6[) are 
most convenient for relating the BS wave function to the four-point Green's function 
and the scattering observables at L — > oo. Closely related observation was obtained 
long time ago by Nishijima, Zimmermann and Hagg who derived the generalized 
reduction formula for local composite fields.^ 

In principle, one may choose any composite operators with the same quantum 
numbers as the nucleon to define the BS wave functional] Different operators give 
different BS wave functions and different AW potentials, although they lead to 
the same observables such as the phase shifts and binding energies. This is quite 
analogous to the situation in quantum mechanics that the unitary transformation of 
the wave function changes the structure of the potential while the observables are not 
modified. A theoretical advantage of our approach based on lattice QCD is that we 
can unambiguously trace the one-to-one correspondence between the AW potential 
and the interpolating operator in QCD. This is in contrast to the phenomenological 
AW potentials where connections to QCD operators are not attainable. 

*' In practice, however, we had better restrict ourselves to consider only local composite opera- 
tors for the nucleon, since it is very difficult, although not entirely impossible, to derive the reduction 
formula for non-local composite operators without violating the causality of relativistic theories. 
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Table I. Two-nucleon asymptotic states classified by the total isospin /, the total spin (s), the or- 
bital angular momentum (£), and the total angular momentum (J) together with some examples 
in low partial waves. 



I 





1 


s 





1 





1 


£ 


odd 


even 


even 


odd 


J 


I 


£ 




£ 


£ 


l±l 


J = 








'So 




3 Po 


J = 1 


'Pi 




3 Si, 3 Dx 




3 Pr 




J = 2 




3 D 2 








3 P 2 , 3 F 2 


J = 3 






3 D 3 , 3 G 3 




3 F a 




J = 4 




3 G 4 




X G 4 




3 F 4 , 3 H 4 

















§4. General form of the NN potential 

In the previous section, we illustrated the procedure to define the potential 
between the neutron and the proton, which has spinor indices a, f3, 7, 5 running from 
1 to 4. In order to derive the general structure of the NN potential at low energies, 
we restrict ourselves to consider only the upper components of these spinor indices 
in the following sections. 

4.1. Symmetry of the two nucleon system 

It is useful to classify the asymptotic two-particle states by the orbital angular 
momentum (£), the total spin (s) and the total angular momentum (J) together with 
the total isospin /. Using the standard notation, 2s+1 £j, and taking into account 
constraints due to Pauli principle, we have the well-known relations given in Table 
I. 

4.2. Okubo-Marshak decomposition 

The general form of the NN potential in the two-component spinor space has 
been classified by Okubo and MarshakPS We leave the derivation in Appendix 
B and recapitulate only the results here. By using the helmiticity, translational 
invariance in space and time, Galilei invariance, rotational invariance, parity and 
time-reversal invariance, fermi statistics and isospin invariance, the potential has a 
general decomposition 

V = Y / V I (r,v,a l ,a 2 )PJ, (4-1) 
1 

V 1 = Vi + Vl (a, ■ <r 2 ) + ±{V$, S 12 } + v[ s L-S+ i{y£, P 12 } + \{V&, W 12 }, 

(4-2) 



10 S. Aoki, T. Hatsuda and N. Ishii 

where P/ is the projection operator to the iso-singlet (/ = 0) and iso-triplet (J = 1): 

P5 = \-ri-r 2 , Pl = ^ +T1 . T2 . (4-3) 

Also, we define 

S12 = 3(<xi • f)(er 2 ■ f) - <ri ■ er 2 , (4-4) 

S = -(o-i + <r 2 ), L = rxp, (4-5) 

Fi2 = (o-i-v)(o- 2 -t;), (4-6) 

VFi2 = Ql2-^(o-l-cr 2 )i: 2 , (4-7) 

Q 12 = -{<7i-.L,<7 2 -£}, (4-8) 

with v = p/ 1 jj,. The anticommutators in Eq. (|4-2j) are necessary to make the po- 
tential hermitian, since Si 2 , P± 2 , W\ 2 do not commute with the scalar potentials 
Vi(r 2 ,v 2 ,L 2 ) (A = Q,a,T,LS,P,W). 

If we keep the terms only up to the first order in v, we obtain the conventional 
form of the potential at low energies commonly used in nuclear physics: 

V 1 = Vi{r) + V» (<Ti • <x 2 ) + Vi{r) S 12 + V[ s {r) L-S + 0(v 2 ), (4-9) 

or in a more conventional notation, 

V = V c (r) + V T (r)S l2 + V LS (r)L ■ S + 0(v 2 ), (4-10) 
= Vb(r) + V a {r)(ai ■ <t 2 ) + V T {r){n ■ r 2 ) + V r7T (r){a 1 ■ a 2 ){r 1 ■ r 2 ) 
+ [^T0(r) + ^Tr(r)(Ti-T 2 )] S 12 

+ [V L so(r) + Vbsr(r)(Ti • r 2 )] L ■ S + 0(v 2 ). (4-11) 

The central and tensor potentials, and Vy, in Eq. (|4-10|) are the leading-order 
(LO) terms of O(v ) in the velocity expansion, while the spin-orbit potential, Vls 
is the next-to-leading-order (NLO) term of 0{v). 

4.3. Determination of the NN potentials 

For given /, s and J, the matrix elements of the LO and NLO potentials up 
to 0(v) in Eq. (|4-11|) have the following structure (see Appendix C and also see 
Ref-EHD): 

V^r; 1 J J ) = V I (r) + V*(r), (4-12) 
V\r- %) = Vi(r) - W*{r) + 2^(r) - V&(r), (4-13) 
v // r . * U s _ ( Vl_(r) VL + (r) \ 

V (r, (J T l)j)- ^ F /_ (r) y / +(r) J, (4-14) 



with 



KL(r) = ttf (r) - 3F» - ^— ]V/(r) + ( J - l)V&(r), (4-15) 
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Vl + (r) = V I (r)-3Vi(r) 



2(J + 2) 



Vl+(r) = Vl_{r) = 6 



2J+ 1 



y/(r)-(J + 2)y/ s (r) 



2J + 1 

rl=0,l 



l#(r). 



(4-16) 
(4-17) 



There are 8 unknown functions, V ~£ ST , while we have 4 (2) diagonal and 1(0) 
off-diagonal matrix elements at each J for J > (J = 0) as seen from Table I. On 
the lattice, it is relatively unambiguous to extract information for i = 0,1,2,3 = 
S,P,D,F using the irreducible representations of the cubic groupP^ 1 Then, at most 
16 independent (14 diagonal and 2 off-diagonal) information as seen in Table I are 
obtained for 8 unknowns Vj[(r), so that each Vj[(r) can be determined in two different 
ways. 

4.4. Long range part of the potential 

In QCD with dynamical quarks, the lightest hadron is the pion. Therefore, the 
longest range interaction between the nucleons is dictated by the one-pion-exchange 
potential (OPEP). For later purpose, let us here summarize several features of OPEP 
with special care about its chiral behavior. 

First of all, the equivalence theorem implies that the pseudo-scalar irN coupling 
g n N(— 14.0) and the pseudo- vector coupling f w N at low energy are related through 
fnN = 2% • This is simply obtained by kinematics. On the other hand, chiral 
symmetry leads to the Goldberger-Treiman (GT) relation, ~ jr-, where g A {— 
1.27) is the nucleon axial-charge and F n (~ 93 MeV) is the pion decay constant. 

With these relations, the OPEP reads 



VoPEp(r) 



ft 



N 



4tt 
,2 



Tl • T 2 )(<Ti • Vl)(<T 2 • v 2 )- 
2 



4vr \2M N ) 3 



(a i ■ a 2 ) + S 12 1 + 



+ 



m n r mir 2 



9 A ( m n \ {Ti ■ r 2 ) 



4vr V 2 R 



chiral limit 



[r\ ■ T 2 ) 



£Ti • <T 2 ) + S12 

Sl2 



3 3 
1 + + 



m n r mlr 2 



r 

-m n r 



(4-18) 
(4-19) 

(4-20) 
(4-21) 



Here we have used the equivalence theorem to obtain Eq. (I4-19P from Eq. (14-180 and 
use the GT relation to obtain Eq. (|4~20D from Eq. gd9j). g A and F n in Eq. <^TDi 
are the values in the chiral limit. 

In quenched QCD without dynamical quarks, there arises a dipole ghost in the 
flavor-singlet channel (the 77-channel in the case of two flavors) which couples to the 
nucleonsPJ'ES The 7/-propagator in the quenched approximation is written as 

iM 2 (q) 



D v (q) 



+ 



q 2 — m 2 + ie (q 2 — m 2 + ie) 2 



(4-22) 



where M^q) 



a>oq 2 with mo and «o being ghost parameters. The second 
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term is the dipole ghost corresponding to the hairpin diagram with quark-line dis- 
connected. Then the NN potential from the r\ exchange reads 3 ^* 



V v (r) 



!^.V)(<x 2 .V) 

-2 



;i-a ) + M 2 (m 7r ; 



d 



o 



rat 



(4-23) 



9 v n ( m n \ (1 - ap) 
4vr \2M N ) 3 

9 2 r,N ( m n \ 2 / M 2 (m^) \ 1 



3 3 
Ol • cr 2 ) + 5i2 ( 1 + + 



m T r mlr 2 



;cr.l. ■ cr 2 ) | 1 — ) + S12 ( 1 + 



(4-24) 



where f^N (gr]N) is the pseudo-vector (pseudo-scalar) coupling of the flavor-singlet r] 
to the nucleon. Its magnitude does not necessarily be as large as the ttN coupling.® 
Note that the long range part of the potential has exponential fall-off instead of the 
Yukawa- type because of the dipole-term in Eq. ()4-22p . 

Let us define a ratio IZ13 between the central potential in the spin-singlet channel 
and that in the spin-triplet channel, 



n 



13 



Vc{r\ 'Sp) 
V c (r; 3 Si) 



+1 (one-pion-exchange) , 
— 3 (one-ghost-exchange) . 



(4-25) 



Since we have (eri • <T 2 ) S pi n _ sin gi c t = -3, (crx ■ 0-2) spin-triplet = +1 and the similar 
relations for the isospin, the large r behavior of IZ13 has different sign and magnitude 
between the one-ghost-exchange and one-pion-exchange. Therefore 7^13 can be used 
as a tool to identify the ghost contribution at large distance as will be discussed in 



§5. Central and tensor forces in lattice QCD 



5.1. BS wave function on the lattice 

To define the BS wave function on the lattice with the lattice spacing a and the 
spatial lattice volume L 3 , we start from the four-point correlator, 



Qap(x,y,t- t ; J 



(0 \n p (y, t) Pa (x,t)J pn (t ; J p )\0) (5-1) 

00 

^2A n {0\n p {y)p a (x)\E n ) e -^(*-*>), (5-2) 

n=0 

> A)Vv(r;J P )e-^-'«), (5-3) 



with the matrix element A n = (E n \J ' pn (0)|0) . The states created by the source 
J pn have the conserved quantum numbers, (J, J z ) (total angular momentum and its 
z-component) and P (parity). For studying the nuclear force in the J p = + ( 1 So) 
channel and the J p = 1 + ( 3 Si and 3 Di) channel, we adopt a wall source located at 
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t = to with the Coulomb gauge fixing only at t = to: 



j pn (t o; j p ) = [ P r (^r^j . ( 5 - 4 ) 

where Pa all (to) and n^ all (to) are obtained by replacing the local quark fields q(x) and 
g(y) in Eqs. (|3-5p and (j3-6p by the wall quark fields, 

^(to) = to). (5-5) 



By construction, the source operator Eq. (!5-4p has zero orbital angular momentum 
at t = to, so that states with fixed (J,J Z ) are obtained by the spin projection with 
(s, s z ) = (J, J z ), e.g. P^~ 0) = (0-2)^ and P^~ 1,S "~ Q ' = (<ri)p a . Note that the £ and 
s are not separately conserved: Therefore, the state created by the source J pn (to; 1 + ) 
becomes a mixture of the i = and t = 2 at later time t. 

The BS wave function in the orbital S-state is then defined with the projection 
operator for the orbital angular momentum (PW) and that for the spin (P( s )): 



V^So) = P^=°)p( s = V(r;0+) = ^E P £ =0) ^(^;0+), (5-6) 

<?eO 

V(r; 3 S X ) = P^=°)p( s=1 ^(r; 1+) = 1 £ P£ =1 W<rV 1+). (5-7) 



24 

9 eO 



Here the summation over g E O is ta ken for the cubic transformation group with 24 
elements to project out the S-stateEl III 



5.2. Asymptotic momentum 

The asymptotic momentum A; for the S-states is obtained by fitting the BS wave 
function ip{r) with the Green's function in a finite and periodic boxPS 1 

1 x - i{2ir/L)n-r 

which satisfies (V 2 + k 2 )G(r;k 2 ) = — 5\ at (r) with <5i a t(r) being the periodic delta- 
function. In the actual calculation, Eq. (15-80 is rewritten in terms of the heat kernel 
K, satisfying the heat equation, dtK,{t,r) = V 2 /C(i,r) with the initial condition, 
1C(t — > + ,r) = 5i a t(?") (see Appendix ID1 for the detail). The fits are performed 
outside the range of the NN interaction determined by V 2 ip(r)/ip(r)W^ 



More precisely, this projection picks up an Af state, which contains not only an £ = 
component but also the higher orbital waves with I > 4. Latter contributions, however, are expected 
to be negligible at low energy. 

**' Note that p^ = °)p( s=0 > in Eq. (|5-6[) is a redundant operation, since we have already prepared 
J p = + state by the wall source J pn (to;Q + ) which allows only the 1 So channel. Also, p( s=1 ) in 
Eq. (|5-7p is a redundant operation, since the J p = 1 + state prepared by the wall source J pn (to', 1 + ) 
allows only the spin-triplet state. 
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5.3. Effective central potential at low energies 

In the S-states at low energies, the effect of the velocity dependent terms in 
Eq. (14- lip is supposed to be small compared to the velocity independent terms, so 
that it is convenient to define the "effective" central potential V^ ff (r)p2J 

^M=* + -^. (5-9) 
mN ip{r) 

As long as we keep only the LO terms of the velocity expansion in Eq. (|4-10p . 
V^, s (r; 1 So) is equivalent to Voir; 1 So), while Vq (r; 3 Si) differs from Vc{r; 3 Si) due 
to the higher order effects from the tensor potential. One can also study the validity 
of velocity expansion in Eq. (|4-10p by calculating VQ S (r; l So) for different energies 
E (see Ml. 



5.4. Scattering lengths 

The NN scattering lengths for the S-states can be deduced from Liischer's for- 
mu lajI3J33 

k cot 5 (k) = -^Z 00 (l;q 2 ) = - + 0(k 2 ), (5-10) 

where Zoo(l;o 2 ) with q = |— is obtained by the analytic continuation of the gen- 
eralized zeta-function Zoo(s;q 2 ) = —r= ^™ez 3 ( n2 ~ Q 2 )~ s defined for Re s > 3/2. 



(See also Ref. 38) for more general considerations.) In this formula, the sign of the 
S-wave scattering length ao is defined to be positive for weak attraction. 

5.5. Decomposition into central and tensor potentials 

Although the tensor force at long distance is dominated by the one-pion ex- 
change, its spatial structure at medium and short distances is not well understood 
theoretically nor well determined phenomenologically. Therefore, it is quite impor- 
tant to extract it from lattice QCD. 

In the LO of the velocity expansion in Eq. ()4-10p . only the central potential Vc(r) 
and the tensor potential Vt{t) are relevant: The central potential acts separately 
on the S and D components, while the tensor potential provides a coupling between 
these two. Therefore, we consider a coupled-channel Schrodinger equation in the 
J p = 1 + channel, 39 ' in which the BS wave function has both S-wave and D-wave 
components: 

[H + V c {r) + V T {r)S l2 )^{r- 1+) = E^(r; 1+). (5-11) 

The projections to the S-wave and D-wave components similar to Eq. ()5-7|> read 

V^ap = P {i=0) ipa(i(r; 1+), (5-12) 
Qifrap = (1 - P (e=0) )i>ap(r; 1 + ). (5-13) 

Note that both Vijjap and Qtpa/3 contain additional components with £ > 4 but they 
are expected to be small at low energies. 
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By multiplying V and Q to Eq. (|5-lip from the left and using the fact that Hq, 
Vc(r) and Vr(r) commute with V and Q, Eq. (|5-11|) splits into two equations, 

H [Pi/>](r) + V c (r)[PiJ>](r) + V T {r)[V S 12 ip]{r) = E[V^](r), (5-14) 
H [Qip](r) + V c (r)[Qi>](r) + V T (r)[QS 12 ip](r) = E[Q^]{r), (5-15) 

where we have suppressed the spin indices, a and /3, for simplicity. 

By picking up (a, f3) = (2, 1) component of these two equations, we arrive at 

V c (r) = E- _l_([Q5i2^] 2 i (r)H Q [V^) 2 i{r) - [PS 12 iphi(r)H [Qrphi(r)), (5-16) 
Ayr) v ' 

V T (r) = ^y([GV]2i(r)fr [^] 2 i(r) - [PV>]2i(r)#o[Q^2i(r)), (5-17) 

A(r) = [T^}2i(r)[QS 12 ^hi(r) - [Qiphi(r)[PS 12 iP] 21 (r). (5-18) 

§6. Numerical results in quenched QCD 

6.1. Setup of the lattice simulations 

We employ the standard plaquette gauge action on a 32 4 lattice with the bare 
QCD coupling constant (3 = 6/g 2 = 5.7. The corresponding lattice spacing is deter- 
mined as 1/a = 1.44(2) GeV (a ~ 0.137 fm) from the p meson mass in the chiral 
limit.® The physical size of our lattice then reads L ~ 4.4 fm. As for the fermion 
action, we adopt the standard Wilson quark action with the hopping parameter 
(« = 0.1640,0.1665 and 0.1678), which controls the quark masses. The periodic 
boundary condition is imposed on the quark fields along the spatial direction, while 
the Dirichlet boundary condition is imposed along the temporal direction on the 
time-slice t = 0. The wall source is placed on the time-slice at to /a = 5 after the 
Coulomb gauge fixing at t = to. 

To generate the quenched gauge configurations, we adopt the heatbath algorithm 
and sample configurations are taken in every 200 sweeps after skipping 3000 sweeps 
for thermalization. The number of sampled gauge configurations -/V con f, the pion 
mass m n , the rho- meson mass m p and the nucleon mass are summarized in 
Table HH For k = 0.1678, we have removed 28 exceptional gauge configurations from 
the sample. 

The BS wave functions are measured at (t — to)/a = 7, 6, 5 for k = 0.1640, 0.1665, 
0.1678, respectively. These values of t — to are determined by studying the ground 
state saturation in the iViV potentials as discussed below. We employ the nearest 
neighbor representation of the discretized Laplacian as X7 2 f(x) = Ylf=i{f( x + an i) 
+ f{x — arij)} — 6f(x), where rii denotes the unit vector along the i-th coordinate 
axis. BS wave functions are fully measured for r < 0.7 fm, where rapid change of 
the AW potential is expected. Since the change is rather modest for r > 0.7 fm, the 
measurement of BS wave functions has been restricted on the coordinate axes and 
their nearest neighbors to reduce the computational cost. 
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Table II. Summary of the hopping parameter re, the pion mass m^, the rho-meson mass m p , the 
nucleon mass jtin, the time-slice (t — to)/a at which BS wave functions are extracted, the spatial- 
slice R/a above which the NN potentials are inactive, and the number of gauge configurations 
A'conf with exceptional configurations being removed. The lattice spacing is a ~ 0.137 fm. Some 
numbers are updated from Tables 1 and 2 of Ref. I21|l . 



re 


[MeV] 


m p [MeV] 


mjv [MeV] 


(t - t ) /a 


R/a 


Neon! 


0.1640 


731.1(4) 


990.3(13) 


1558.4(63) 


7 


11 


1000 


0.1665 


529.0(4) 


894.3(28) 


1333.8(82) 


6 


11 


2000 


0.1678 


379.7(9) 


837.9(21) 


1196.6(83) 


5 


12 


2021 




0.0 0.5 1.0 1.5 2.0 
r[fm] 

Fig. 1. The NN wave functions in 1 S and 3 Si channels for m n = 529 MeV (re = 0.1665). The 
inset is a three-dimensional plot of the wave function ip(x, y, z = 0; 1 So). 



6.2. BS wave functions in the S- state 

Figure [1] shows the BS wave functions in 1 So and 3 Si channels for k = 0.1665. 
The wave functions are normalized to be 1 at the largest spatial point r = 2.192 fm. 

Figure[2ja,b) show the fitting of the wave function in the interval R/a < r/a < 16 
using Eq. (|5-8|) . This leads to the values of the effective energy E = k 2 /m^ in Table 
im The value of R is determined from the ground state saturation of the potential 
as discussed below. 

6.3. Effective central potential 

Shown in Fig. [3] are the reconstructed effective central potentials in the 1 So and 
3 Si channels for k = 0.1665 with the formula Eq. ()5-9j) . The overall structures of the 
potentials are similar to the known phenomenological NN potentials discussed in 
£JH namely the repulsive core at short distance surrounded by the attractive well at 
medium and long distances. From this figure, we find that the interaction between 
the nucleons is well switched off for r > 1.5 fm, so that we chose R/a = 11 (for 
m n = 731, 529 MeV) and R/a = 12 (for m T = 380 MeV) as given in Table IE In 
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0.0 0.5 1.0 1.5 2.0 

r[fm] 

Fig. 2. (a) The fit of the NN wave functions for m w — 529 MeV in the 1 So channel using the 
Green's function in the fit range 11 < r/a < 16. (b) Similar fit for the NN wave functions in 
the 3 Si channel. 



both cases, the condition R < L/2 = 2.2 fm is satisfied. 

To check the stability of these potentials against the time-slice adopted to define 
the BS wave functions, we plot the i-dependence of the 1 So potential for several 
different values of r as shown in Fig. U for m n = 529 MeV: In this case, choosing 
(t — to) /a = 6 to extract Vc(r) would be good enough to assure the stability within 
the statistical errors. The time-slices chosen for other cases by the same procedure 
are given in Table HT1 

6.4. Quark mass dependence of the central potential 

In Fig. [5l we compare the NN central potentials in the 1 So channel for three 
different quark masses. As the quark mass decreases, the repulsive core at short 
distance and the attractive well at medium distance are enhanced simultaneously. 
This feature can be also seen in Fig. E^a) where r 2 Vc(r), which appears in the 
quantum mechanical matrix elements, is plotted. To study the relative magnitude 
of the repulsion and the attraction, we define the following volume integrals of the 
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Fig. 3. The effective central potentials in the 1 So channel and in the 3 Si channel for ro, = 529 
MeV. 



potential and plot them in Fig. E^b): 

pro m 
Jl = / r 2 V c (r)dr, I 2 = r 2 V c {r)dr. (6-1) 

JO Jr 

Here ro (~ 0.5 fm) is the first nodal point where r 2 Vc(r) changes sign from positive 
to negative, and ri is the point at which r 2 Vc(r) becomes essentially zero within the 
statistical errors. The error bars in Fig. E^b) reflect the uncertainties of ro,i as well 
as those from the spline curve fit of the data. The comparison of 1%, I 2 and I\ + I2 
implies that (i) both repulsion and attraction increase in magnitude as quark mass 
decreases, and (ii) there is a large cancellation between the repulsion and attraction, 
and (iii) there is a net attraction increasing as the quark mass decreases. 

6.5. Dipole ghost in the central potential 

To check if there is an evidence of the exponential tail from the dipole ghost in the 
long range part of the effective central potentials, the ratio IZ13 given by Eq. (j4-25j) 
is plotted in Fig. \7\ as a function of r for the lightest quark mass, = 380 MeV. 
Within the statistical errors, there is no sign that 7^13 — > — 3 for r > 1 fm, so 
that possible ghost contamination is small in our results with relatively heavy quark 
masses. The figure also shows that 7^13 is rather close to +1 for r > 0.7 fm. This 
does not necessary implies that the OPEP is seen: as long as there are spin-isospin 
independent attraction such as originating from the two-pion-exchange potential, it 
also leads to 7^13 ~ 1. 
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Fig. 4. The t-dependence of the potential at r = 0,0.14,0.19, 1.37, 2.19,0.69 fm from top to bottom 
for the 1 So channel at = 529 MeV. 



Table III. Effective center of mass energies E = fc 2 /mjv obtained from the asymptotic momenta for 
different quark mass, ao's are the associated scattering lengths obtained from Liischer's formula 
Eq. (JSIIOl). 



m n [MeV] 


E^Sa) [MeV] 


JS( 3 Si) [MeV] 


ooCSo) [fm] 


ao( 3 Si) [fm] 


731.1(4) 


-0.400(83) 


-0.480(97) 


0.115(26) 


0.141(31) 


529.0(4) 


-0.509(94) 


-0.560(114) 


0.126(25) 


0.140(31) 


379.7(9) 


-0.675(264) 


-0.968(374) 


0.153(66) 


0.230(101) 



6.6. NN scattering lengths 

As we found in Fig. the central potential multiplied by r 2 shows a net attrac- 
tion as a result of the large cancellation between the short range repulsion and the 
medium range attraction. This attractive nature of the potential can be quantified 
by the scattering length ao defined from Liischer's formula, Eq. (|5-10l) . together with 
the asymptotic momentum k obtained from Eq. (|5-8p Fl 

The results of ao are summarized in the last two columns in Table IIIII where 
0(k 2 ) correction on the right-hand side of Eq. (|5- 101) is assumed to be small for 
the present energy E = k 2 /mj^. In Fig. [U the scattering lengths for 1 So and 3 Si 
channels are shown as a function of m 2 . Although there is a small attraction which 
increases as m n decreases in both channels, the absolute magnitudes of ao are much 

*' If the net interaction is small in the infinite volume limit, the volume integral of the 
potential and the scattering length are related in the Born approximation as, a ^' cak - cou P lm g 
~ — mjv J Vc(r)r 2 dr. 
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Fig. 5. The central potentials in the 1 So channel for three different quark masses. 



smaller than the experimental values at the physical point: ( 1 So) ~ 20 fm and 

4 cxp) ( 3 Si) ~ -5 fm at ml = 0.018 GeV 2 . 

The above discrepancy is partly attributed to the heavy quark masses employed 
in our simulations: If we can get closer to the physical quark mass in full QCD 
simulations, there should arise the "unitary region" where the NN scattering length 
becomes singular and changes sign. This was first noted in clear terms by Kura- 
mashP^J anc [ was later elaborated in Refs. H9|) and FIT]) by using chiral perturbation 
theory. The singularity is associated with the formation of the di-nucleon bound 
state, so that the NN scattering length becomes a non-linear function of the quark 
mass in the unitary region. As suggested in Ref . I23j) by using the one-boson-exchange 
model with the quark-mass dependence of the hadron masses taken from the lattice 
QCD data, the size of the unitary region could be narrow, which implies that the 
scattering lengths at the heavy quark masses adopted in our simulation can be as 
small as the values in Fig. El 

Unlike the scattering length, the NN potential would not have singular behavior 
in the unitary region as expected from the well-known quantum mechanical examples 
such as the low-energy scattering between ultracold atoms. Also, the effective range 
parameter would be a rather smooth function of the quark mass. To check these 
points in QCD, it is important to study the AW potential, the scattering length 
and the effective range simultaneously in the full QCD simulations which allow us to 
approach small quark masses without quenched artifact. Studies along this direction 
is now underway^ 1 and will be reported elsewhere. 
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Fig. 6. (a) The central potentials with r 2 multiplied in the 1 So channel for three different quark 
masses, (b) Comparison of the attractive part and repulsive part of the potential in terms of 
the volume integral in the So channel. 



6.7. BS wave function in the ~D- state 

In Fig. [Sfa), we show the 3 Si and 3 Di components of the BS wave functions 
obtained from the J p = 1 + , J z = M = state for m n ~ 529 MeV, according to the 
procedure given in £ )5,5l To reduce the computational cost, the points are restricted 
on the coordinate axes and their nearest neighbors for r > 0.7 fm, whereas all points 
are calculated for r < 0.7 fm. 

Note that the 3 Di wave function as a function of r is multivalued due to its angu- 
lar dependence. Since (a, j3) = (2, 1) spin component of the D-state wave function for 
J p = 1 + , M = is proportional to the spherical harmonics Y2o(@, <j>) oc 3 cos 2 6 — 1, 
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Fig. 7. The ratio of the central potentials defined in Eq. (|4-25[) for the lightest quark mass, 
380 MeV. 




m* 2 [GeV 2 ] 

Fig. 8. Scattering length ao in the 1 So and 3 Si channels for three different quark masses obtained 
in the quenched QCD simulations. 



it is a good consistency test to check if the multivaluedness can be absorbed by this 
angular dependence. Shown in Fig. Mjo) are the same BS wave functions as Fig. E(a) 
with the angular dependence in the D-state assumed to have this spherical harmon- 
ics form. It is clear that the multivaluedness is nicely removed, and thus it is certain 
that we indeed extracted the D-state wave function on the lattice. 

6.8. Tensor force and its quark mass dependence 

Shown in Fig. [10] are the central potential Vc{r) and tensor potential Vt(t) 
together with effective central potential Vq(t) in the 3 Si channel. (As mentioned 
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Fig. 9. (a) (a, ft) = (2, 1) components of the S-state and the D-state BS wave functions projected 
out from a single state with J p = 1 + , M = 0. (b) The same data with the spherical harmonics 
components are removed in the D-state. 



before, we consider only the LO terms of the velocity expansion here by assuming 
that the NLO term (the spin-orbit potential) and higher order terms are negligible 
at this low energy.) 

Note that VQ S (r) contains the effect of Vr(r) implicitly as higher order effects 
through the process such as 3 Si — ?- 3 Di — > 3 Si. In the real world, VQ ff (r) is expected 
to acquire sufficient attraction from the tensor force. This is the reason why bound 
deuteron exists in the 3 Si channel while the bound dineutron does not exist in the 
1 So channel. Now, we see from Fig.[TU]that the difference between Vc(r) and V^ ff (r) 
is still small in our quenched simulations due to relatively large quark masses. This 
is also consistent with the results of the small scattering length shown in Fig. [8] 

The tensor potentials Vr(r) in Fig. [10] are negative for the whole range of r 
within statistical errors and have a minimum at short distance around 0.4 fm. If the 
tensor force receives significant contribution from the one-pion exchange as expected 
from the meson theory, Vr(r) would be rather sensitive to the change of the quark 
mass. As shown in Fig. [TT] it is indeed the case: Attraction of Vr(r) is substantially 
enhanced as the quark mass decreases. A phenomenological fit of the tensor force 
taking into account this physics will be given later. 

As discussed in ^6.5t the ratio IZ13 of the effective central potentials in the 1 So 
and 3 Si channels is close to unity for r > 0.7 fm so that we do not see evidence 
of the dipole ghost (quenched artifact) in the long range part of the potential with 
our relatively heavy quark masses. However, this does not necessary imply that the 
OPEP is seen in the effective central potentials: If the OPEP dominates at long 
distances, Eq. (|4-19p immediately implies that the magnitude of the tensor potential 
is always larger than the central potential at long distances. Since this is not seen 
in Fig. [10] within the statistical errors, it is unlikely to interpret the attraction of 
V(f(r) at 0.5 fm < r < 1 fm as the evidence of OPEP. 

A technical comment is in order here. Since we use the (a,/3) = (2,1) spin 
component of Eq. (j5-16p . the second equation vanishes at r oc (±1,±1,±1). This 
is because the spin (2, 1) component of the D-state wave function is proportional to 
Y2o(0,(f)) oc 3cos 2 # — 1 which vanishes at r oc (±1,±1,±1). Although these points 
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Fig. 10. The central potential Vc{r) and the tensor potential Vr(r) obtained from the J p = l" 1 
BS wave function at = 529 MeV. 
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Fig. 11. Quark mass dependence of tensor force. The lines are the four-parameter fit using the 
one-p-exchange + one-pion-exchange with Gaussian form factors. 



are removed from our plots, statistical error is accumulated in the neighborhood of 
these points. (For instance, see the points at r ~ 0.5 fm in Figs. [TO] and [TTJ) A 
resolution of this problem by combining the data with other spin components will 
be reported in the future publication. 

The central and tensor potentials obtained from lattice QCD are given at dis- 



Lattice Nuclear Force 



25 



crete data points. For practical applications to nuclear physics, it is more useful 
to parametrize the lattice results by known functions. We have tried such a fit 
for Vr(r) under the assumption of the one-/>exchange + one-pion-exchange with 
Gaussian form factors: 



where, ^1,2,3,4 are the fitting parameters while m p (m n ) is taken to be the />meson 
mass (the pion mass) calculated for each quark mass. At this moment, it is hasty to 
extract physical quantities from the fit such as the meson- nucleon coupling constants: 
Nevertheless, it may be worth mentioning that the pion-nucleon coupling constant 
extracted from the parameter 63 in the case of the lightest pion mass (m n = 380 
MeV) reads gl N /{±it) = 12.1 ± 2.7 which is encouragingly close to the empirical 
value. We have tried similar fits for the central potential with the phenomenological 
repulsive core with a Gaussian form and the meson-exchange potential with form- 
factors: The results are still not stable enough due to the statistical errors of the 
lattice data. 

6.9. Velocity dependence of the potential 

So far we have considered the potential determined from the lattice data taken 
almost at zero effective energy E ~ MeV (see Table UTTJ). If the local potential 
determined from the other energies has different spatial structure, it is an indication 
that there are velocity dependent terms as discussed in ^2.1i 

A lattice QCD analysis on the velocity dependence has been recently carried 
out by changing the spatial boundary condition of the quark field from the periodic 
one to the anti-periodic one, so that the effective center of mass energy is increased 
to E ~ 3(-7t/ L) 2 jrnjq ~ 50 MeVP^l 1 The result shows that the central and tensor 
potentials do not show modifications for every r within the statistical errors: Namely, 
the non-locality of the potential with our choice of the interpolating operator is small 
and the potentials shown in the present paper can be used in the energy region at 
least up to E ~ 50 MeV without significant modifications Detailed account of the 
above result is beyond the scope of this paper, and will be reported elsewhere. 



In this paper, we have discussed the basic notion of the nucleon-nucleon potential 
and its field-theoretical derivation from the equal-time Bethe-Salpeter wave function 
in QCD. By construction, the non-local potential defined through the projection of 
the wave function to the interaction region (the inner region) correctly reproduces 

*' An investigation based on integrable models suggests that potentials derived from the BS 
wave functions with local operators in these models are slowly varying functions of energy (veloc- 





§7. Summary and concluding remarks 



ity)P 



26 



S. Aoki, T. Hatsuda and N. Ishii 



the asymptotic form of the wave function in the region beyond the range of the nu- 
clear force (the outer region). Thus the observables such as the phase shifts and the 
binding energies can be calculated after extrapolating the potential to the infinite 
volume limit. Non-locality of the potential can be taken into account successively 
by making its velocity expansion, which introduces the velocity-dependent local po- 
tentials. The leading-order terms of such velocity expansion for the nucleon-nucleon 
interaction are the central and the tensor potentials. 

As an exploratory study to test how this formulation works, we have carried 
out quenched lattice QCD simulations of the two-nucleon system in a spatial box of 
the size (4.4 fm) 3 with the quark masses corresponding to m n = 380, 529, 731 MeV. 
We found that the NN potential calculated on the lattice at low energy shows all 
the characteristic features expected from the empirical NN potentials obtained from 
the experimental NN phase shifts, namely the attractive well at long and medium 
distances and the repulsive core at short distance for the central potential. As for the 
tensor potential obtained from the coupled channel treatment of the 3 Si-state and 
the 3 Di-state in the BS wave functions on the lattice, we found appreciable attraction 
at long and medium distances and a moderate repulsion at short distance. 

As the quark mass decreases, the repulsive core and attractive well in the central 
potential, and the attractive well in the tensor potential tend to be enhanced. Also, 
we found net attraction in both 1 So and 3 Si channels after the cancellation of the 
repulsive core and the attractive well. The absolute magnitudes of the scattering 
lengths are still much smaller than the physical values due to the large quark mass 
in our simulation. Phenomenological fit of the tensor potential strongly suggests the 
existence of the one-pion-exchange contribution in its long range part. 

There are a number of directions to be investigated on the basis of our approach 
as listed below: 

1. Determination of the velocity dependence is important in deriving the NN po- 
tentials which can be used for the wide range of scattering energies. Studies 
along this line using the anti-periodic boundary condition in the spatial direc- 
tion has been already started^ as mentioned in §6.91 

2. To derive the realistic NN potentials on the lattice, it is necessary to carry out 
full QCD simulations with dynamical quarks. Studies along this line with the 
use of the (2+l)-flavor QCD configurations with the Wilson fermion generated 
by PACS-CS Collaboration® is currently under wayP^ 1 

3. The hyperon-nucleon (YN) and hyperon-hyperon (YY) potentials are essential 
for understanding the properties of hyper nuclei and the hyperonic matter inside 
the neutron stars. However, the experimental scattering data are very limited 
due to the short life-time of hyperons. On the other hand, the NN, YN and 
YY interactions on the lattice can be treated in the same manner by changing 
only the quark flavors. Recently, the EN potential in quenched QCD 241 and 
the AN potential in quenched and full QCD 4 ^ are examined as a first step 
toward systematic derivation of the hyperon potentials. 

4. The three-nucleon force is thought to play important roles in nuclear structures 
and in the equation of state of high density matter as mentioned in §2.21 Since 
the experimental information is scarce, simulations of the three nucleons on 
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the lattice combined with the method proposed in £ 12.21 may lead to the first 
principle determination of the three-nucleon potential in the near future. 
If it turns out that the program described in this paper indeed works in full 
QCD with realistic quark masses, it would be the promising first step toward the 
understanding of atomic nuclei and neutron stars from the fundamental law of the 
strong interaction, the quantum chromodynamics. 
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Appendix A 

Bethe-Salpeter Wave Function and Its Asymptotic Behaviour 

In this appendix we derive the behaviour of the Bethe-Salpeter (BS) wave func- 
tion at large r, only using the properties of quantum field theories. 

A.l. Unitarity of S -matrix and structure of T -matrix 

We first determine the structure of the AW scattering T-matrix below the pion- 
production threshold. Due to the unitarity of the S'-matrix. S'S = 1 with S = 1+iT, 
we obtain 

(f\T\i) - (f\T^\i) = z^(f\T j \n)(n\T\i). (Al) 

n 

In the case of NN scattering in the center of mass frame such that (k a , s a ) + (k b , s b ) — > 
(k c ,s c ) + (k d ,s d ) whe re k a = (e fc , fc), k b = (e fc , -k) and k c = (e p ,p), k d = (e p , -p) 

with Ek = k 2 + m 2 N and e p = Jp 2 + mfy, we write 

m(Pc, S c ,p d , S d \T\p a , S a ,p b , S b )i n = (27r)V 4) (p a + p b - p c - p d )T(p, S c , S d ] k, S a , S b ) . 

(A2) 

Here Si = ±1/2 is a helicity of each nucleon, and k = \k\ = \p\ in the center of mass 
frame. Below the pion production threshold such that 2^J k 2 + m? N < 2mj\r + rn^, 
the sum over intermediate states n in Eq. (IA-11) can be restricted to the NN states 
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due to energy-momentum conservations as 

d 3 p± d 3 p: 



Ei-)h = E / (4%(2^%: Ipi ' si ' P2 ' S2)(pi ' si ' P2 ' S21 - (A ' 3) 



S1,S2 

This leads to 



T{p, s c , s d ; k, s a , s b ) - T f (p, s c , s d ; k, 

Sat Sb) 

A; f 

= i z2 32ir 2 £ k J d ^ q Sc ' Sd] 9 ' Sl ' S2 - )T ^' Sl ' S2 ' fc ' Sa ' ( A ' 4 ) 

where |q| = k and fi q is the solid angle of vector q. Using the angular momentum 
basis p*» 

T(p, s c , s d ; k, s a , s b ) = 4vr ^ — (s c , s d |T J (fc)|s a , s b ) {D J )\, ti[ (Q p )D J M s {Q k ) , 

J, A I 

(A-5) 

with s = s a — Sf, and s' = s c — s^, we obtain 

T J (fc) - [T J ]t(fc) = i -h-[T J ]\k)T J (k). (A-6) 

Here is considered as a 4 x 4 matrix and the Wigner D-matrix D J is defined by 

D J MX {Q) = e-* Ma d J MX ((3)e + * Xa , (A-7) 

where the solid angle is denoted as dft = sin (3 df3 da and d J MX (fi) is the Wigner 
d-matrix. The normalization of the D-matrix is given by 

dn (D J ){ M (f2)D J M , x (f2) = ^L_5 JJ '<W, (A-8) 

where no summation is taken for A. For the NN scattering, with new helicity basis 
such that l + j, +|) ± | — \, — \) and | + 5, — \) ± | — \, +\)-, T J is decomposed into 
two lxl submatrices and and one 2x2 submatrix as^ 

T e=j,s=o °lx2 \ 

0' T/ =J)S=1 0i x2 . (A-9) 

2 Xl 2X 1 T t=jzfl, 8 =l I 

The unitarity condition then gives 

TjL JiS = f Js , T/ =Jt1jS=1 = 0(k) ( Tj -^ fj° )o-\k), (A-10) 



with 



% s = ™p.J>»W sfrSuik), O(k) = ( COS£ ^ ( S - Sin£j ^ ] ) , (A41) 
k y smej(A;) cosej(k) /' 1 ; 

where 5i s {k) is the scattering phase shift, whereas £j{k) is the mixing angle between 
^ = J±l. They correspond to the standard Blatt-Biedenharn eigenphase and mixing 
angle. 
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A. 2. BS amplitude and half off-shell T-matrix 

Let us now consider the Bethe-Salpeter (BS) amplitude for the proton and the 
neutron, defined by 

Vapix, y) = (0\T{np(y) Pa (x)}\p(q, s)n(q', s')) in , (A-12) 

where T represents the time-ordered product. The spatial momentum and the he- 
licity for the incoming proton and those for the neutron are denoted by (q, s) and 
(q',s'), respectively. The single nucleon state is normalized covariantly, 
(Bi(q,s)\Bj(q',s')) = 2e q (27r) 3 5ij5 ss >5 3 (q - q') where B\ = p (proton) and B 2 = 
n (neutron). 

The fields, np(y) and p a (x), are the local composite operators for the neutron 
and the proton whose explicit forms are irrelevant for the following derivation. One 
of the advantages to use local operators is that the standard reduction formula can 
be generalized without much modification as shown by Nishijima, Zimmermann and 
Haag (NZH)P^ In particular, one can define in and out composite fields, n in ( out ) (x) 
and p- m ( ut)( x )i i n a similar way as the elementary field through the Yang-Feldman 
equation as^ 1 

V^Vin(out)(z;) = N(x) - / S , ret(adv )(x - x';m)J(x')d 4 x, 



(A-13) 

where A?" takes either n or p, S^ady) denotes the retarded (advanced) Green's func- 
tion in the free space with the mass m = tun, and the "source" is J(x) = (i$ x — 
m)N(x). The wave function renormalization constant Z is defined as ^/~Zu a {p, s) = 
(0| N a (0)\B(p, s)), where we have the following normalization of the Dirac spinors: 



^2ul(p,s)u a (p,s') = ^2vl(p,s)v a (p,s') = 2s p 5 ss i, 

a a 

u a (p, s)up(p, s) = (rf + m) a p, ^2 S )^(P> s ) 

s s 

Then the NZH reduction formula is summarized as 

\t(0)bI( P ,s) - (-)\°\bUp,s)T(0)\ 

d*xe~ ipx T{OA>(x)}[-^- 1 (p)n(p,s)], 

B out (p,s)T(0) - (-)^T(0)B in (p,s) 

d 4 xe ipx [-iu(p,s)S- 1 (p)]T{N(x)0}. 



(A-14) 
m) a p. (A-15) 



(A-16) 



(A-17) 

Here O is an arbitrary product of operators with the number of fermionic operators 
denoted by \0\, and S~ 1 (p) = — m + iS) is the inverse of the free nucleon prop- 
agator. The asymptotic baryon and anti-baryon operators, -B as (p, s) and D as (p,s) 
(as = in, out) are defined by the Fourier decomposition of N^x), 

d 3 p 



aw*) = E 



(2vr) 3 2e t 



B 3S {p,s)u(jp y s) + ^ Dl{p,s)v{p,s) , (A-18) 
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where the flavor and spinor indices are suppressed. The operator £> as thus defined sat- 
isfies the covariant commutation relation, {B as (p, s), -E>| s (p', s')} = 2e p (2Tr) 3 5 ss >5 3 (p— 
p'), and asymptotic states are defined by \B(p : s)) as = £?l s (p, s)|0). 

By using the NZH reduction formula, we can evaluate our BS amplitude Eq. (|A-12p 

as 

&12(X1,X 2 ) = 

Z ' X jll{^Si e ' i9iXi } Gi2MQuQ2;q 3 ,q4)[-iS- 1 (q 3 )u(3M-iS- 1 (q i M4)U, 

(A49) 

where the four-point Green's function is defined by 

r 4 

Gi 2 Mqi,q2;q3,qi) = / [J {d^e****} (0\^{n2(x2)pi(x 1 )p 3 (x 3 )n 4 (x i )}\0). 



i=i 



(A-20) 



Here, to simplify the notation, we abbreviate the Lorentz indices by the lower- 
case suffixes (1, • • • ,4) with the repeated suffixes being contracted and the state 
labels are abbreviated as the numbers in the parenthesis, e.g. u a (q, s) — > u 3 (3) and 
up{q',s') — > ^(4). The four-point function can be decomposed into the free part 

and the connected part as Gi2 ; 34 = Z 2 (G^I- 3 4 + ^^.34)- The free part reads 

G% A = (2n) s 5 4 ( qi - q 3 )5 4 (q 2 - q^iSiq^iS^U, (A-21) 
whereas the connected part is rewritten with the proper vertex r as 

^12:34 (91, 92; 93, qi) 



(2vr) 4 <5 4 (K-Q) [iS{q 1 )\ lv [iS{q2)]22> H)A'2';3'4'0; slQ) [iS(q 3 )] 3 , 3 [iS(q 4 )] 



4/4. 



(A-22) 

Here we have introduced relative and center-of-mass (cm.) 4-momenta by 

K = q 1 +q 2 , k = (q 1 -q 2 )/2, Q = q 3 + q 4 , q = (q 3 - g 4 )/2. (A-23) 
Then, the i-T-integration in Eq. (|A- 19|) can be carried out to obtain 



*i2(xi,x 2 ) = [^(r) + ^g(r)J e~^ H , (A-24) 

4° 2 ) (r) = Z Ul (3)n 2 (4)e-^, (A-25) 
d 4 k 



4kr) = iZ I e- lkr lS( qi )}n'[S(q2)}22>r V 2>Mk;q\Q)M3)M4), 



(A-26) 

where r = x\ — X2 and R = {x\ + x 2 )/2 are relative and cm. 4-dimensional coordi- 
nates, respectively. Covariant Nambu-Bethe-Salpeter type differential equation can 
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be obtained by multiplying 5 _1 (i^i)5 _1 (^ 2 ) to Eqs. (|A~24"1) - (lA~26l) from the left: 

= iZ J ^e- ikr e^ R r a p. l6 (k; q\Q)u 7 (q, s)u s (q', s'). (A-27) 

In our applications of the AW scattering at low energies, it is useful to consider 
the equal-time BS amplitude (which we call the BS wave function in the text) and 
associated Lippmann-Schwinger type integral equation or the Schrddinger type dif- 
ferential equation. For this purpose, we first carry out the integration over k° in 
Eq. (|A-26p using the explicit form of the free propagator: 



S(p) < 1 ^ 1 



171 + iS I a/3 2£ P 



p° — e p + iS p° + e p — id 

(A-28) 

Since we are interested in the asymptotic form of the wave function at \r\ — > oo below 
pion production threshold, we can pick up only the nucleon pole from S(p) in the 
/c°-integral of Eq. (|A-26|) without loss of generality. Possible poles from r associated 
with the resonance production and with the deuteron bound state, as well as anti- 
nucleon poles in S(p) in Eq. (IA-28j) . modify only the short-distant part of the wave 
function. This does not at all imply that those contributions are not important. 
They do affect the actual values of the phase shifts and mixing parameters and are 
fully taken into account in the definition of our potential, Eqs. (|3-3j) and (|3-4p . 

Using the residue theorem and taking the equal-time limit (xq = yo = t) in the 
rest frame of the two-particles {Q = 0), we end up with the Lippmann-Schwinger 
type equation; 

«Mr,t) = ip a p(r-q,s,s') e~ 2ie "\ (A-29) 
i; { $(r;q,s,s') = Zu a (q, s)up(-q, s')e^ r , (A-30) 

^ a f){r; q, s, s') = ^(r; q, s, s') 

r (o) , e q + e k Tss';sAk;q) 

+ (2vr)3^ (r ' fc ' S ' S) k 2 -q*- l 5 + { ] ' { * > 

Here I(r) originates from the contributions other than the nucleon pole and is an 
exponentially localized function in r below inelastic threshold.'^-' In Eq. (IA-3ip . we 
have defined the half off-shell T-matrix, 

*Ti2;34(fc; q) = «i(l)u2(2)(-t)A2 5 34(A;; g|Q)n 3 (3)n 4 (4), (A-32) 

where the outgoing energy 2e^ = 1\J k 2 + m 2 is not necessary equal to the incoming 
energy 2e q = 2y / q 2 + m 2 . The Schrddinger type differential equation is obtained 
from Eq. (jA-31[) by multiplying q 2 + V 2 , 

(q 2 + V 2 )i/j al3 {r;q,s,s') = - ^ J (^3^S^ r; fc ' g ' S ') - 1 ^f^ ' J Ss';s S '{k; q) + £(r), 

(A-33) 
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with JC(r) = (q 2 + V 2 )T(r). Since the plain wave part of ip a p{r; q, s, s') is projected 
out by the operator (q 2 + V 2 ), the right-hand side of Eq. (|A-33p is exponentially 
localized in r and vanishes for r > R. 

A. 3. Asymptotic BS wave function and the phase shift 

Let us further consider the asymptotic behaviour of ip a p(r;q,s,s') at large r 
to relate it to the scattering parameters (phase shifts and mixing angles) defined 
in §A.1I The derivation of this subsection has been essentially given by Ishizuka in 
Ref.EED- 

To perform the k integration in Eq. (jA 33|) . we introduce the following helicity 
decomposition of the half-off shell T-matrix, 

U + 1 

T SlS2 ;s 3 s 4 (k, q) = 4tt Y, -^—(si,s 2 \T J (k; q)\s 3 , Si ){D J )\ M {Q k )D J Ms ,^ q ), 
J,M 

(A-34) 

u a (k, Sl )u (-k,s 2 )e ik - r = ^^ s (a)f/ Qa (V)^(-V)^[ /siS2;a6 (r, k), 

JM 

(A-35) 

where s = s\ — s 2 , s' = S3 — S4, k = \k\ and q = \q\. The reduced wave function </>^1 
in the 2x2 spinor space labeled by the indices a, b is defined as 



'JM\ 



iA 2 M) = YWMis(r,k)(JM£s\JM\ l \ 2 ), (A-36) 



YjMlsi r ^)=u{kr)Yj M {n r ), Yj s M (n r ) = Y J Y Uz {nr)x(s,s z ){lsl z s z \JM). 

(A-37) 

Note that U aa (V) and Upb(— V) in Eq. (|A-35|) are the 4x2 matrices acting on the 2x2 
matrix so that the Dirac structure u a (k, si)up(— k, s 2 ) is correctly reproduced: 
Explicitly, U(V) = \J (sk + mN)(I 2x2 , —icr • V /{Ek + m N))- Alternatively, one may 
use the Lorentz transformation, u(p, s) = A(p)u(Q, s) to define the reduced wave 
function.® 

Note that ( JM£s\ JMXi\ 2 ) in Eq. (IA-36|) is a transformation function between 
the helicity basis and the orbital-spin basis at fixed J, MPS' Also, Xab( s J s z) in 
Eq. (|A-37jl is a 2 x 2 matrix in the spinor space with total spin s = 1 or and its 
z-component s z , and ji(x) is a spherical Bessel function. Using Eq. (|A-8|) . Eq. (|A-3ip 
for large r becomes 



iP(r;q,s,s') = zY,D J MX (n q )U(V)U(-V)^ JM sAr;q), A = s - s' , (A-38) 

JM 

,M , \ sr^ f°° k 2 dk ,ui . , , e + Sk (ss\T J (k: q)\ss') 

s,s' ft 

(A-39) 
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To evaluate the integral in Eq. (|A-39|) . we use the following formula^'® valid 
for r > R in which f(k)k^ i jo(kr)k 2 dk = is satisfiedPl 

00 k 2 dk j e (kr)f{k) . q (+) 

^-^-ir 1 ^ (?r)/(g) - (A40) 

Here hg\x)(=. je(x)±ing(x)) is the spherical Hankel function with jo(x) = (sinx)/x 

and no(x) = —(cosx)/x, so that hg\qr) represents the spherical out going- wave. 
Then, we obtain 



i>JMss'(.r;q) & s A r ^) + { J2 J^^JMrA r ^ q)(ss\T J (q;q)\ss'} 

s,s' ' 

= E [4isA r ^) A fs';sAQ) ~ <>j»,Ar-<l)l^:^<n , (A41) 



s,s' 

A J {q) = l + i-V—T J ( q ;q), B J ( q ) = -^—T J {q-q), (A-42) 

where ' {r, q) is obtained from 0j]v/ss'( r ' ^) by the replacement ji(kr) — > 

ni(qr),h ( e + \qr). 

Using the explicit form of the T-matrix given in (|A-lljl . we finally obtain 

X-/ =J , s = Xjs, Xi =jTltS=1 = 0(q) ( Xj Q 1 ' 1 xj +11 )°' 1{ql (A ' 43) 
with X being either A or B, and 

ifo(g) = e^ (,?) cosfe(g), B £s (q) = e lS ^ sinfe(g), (A-44) 

= 1 (A . 45) 
B ls (q) tanfe(g) 

We now have shown that the BS wave function in QCD has an asymptotic form 
of the scattering wave of the quantum mechanics at large r. To derive this we have 
only use the unitary of the S'-matrix below the inelastic threshold, and have identified 
the phase of the S'-matrix as the scattering phase shift of the asymptotic BS wave 
function. This observation leads to the important conclusion that the potential 
defined through the BS wave function, by construction, gives the correct scattering 
phase shift at asymptotically large r. 



Appendix B 

Okubo-Marshak Decomposition 



In this appendix, we derive the general form of the NN potential in the space 
of two-component spinors, following the argument by Okubo and Marshak. 25 ' The 

*' Since the nucleons are non-interacting in the asymptotic region, the right hand side 
of Eq. (|A-33[1 is exponentially small for r > R. This gives a little weaker condition that 
So° f(k)ji{kr)k 2 dk = 0, which, together with some properties of the T-matrix, leads to the stronger 
condition used here. 
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general form of the 2-body potential with derivatives reads 

V(n,r2,Vl,V2,(Tl,Cr2,Tl,T2,t), (B-l) 

where v lj2 = p lj2 /m N . 

There are several conditions to be satisfied by V. 

1. Probability conservation: This leads to the hermiticity of the potential: V' = V. 

2. Energy- momentum conservation: The energy conservation demands that the 
potential does not depend on time explicitly. The momentum conservation 
leads the translational invariance of the potential. Thus we have 

V = V(r,v 1 ,v 2 , <ri,<r2, ti,t 2 ), (B-2) 

where r = n — r 2 . 

3. Galilei invariance: The potential is assumed to be independent of the center of 
mass momentum of the two-body system, which leads to 

V = V(r,v,tr 1 ,tr 2 ,T 1 ,T 2 ), (B-3) 

where v = p/fi = (p 1 -p 2 )/(2[i) = v 1 - v 2 . 

4. Conservation of total-angular momentum: The total angular momentum is de- 
fined as J = S + L with 

S = -(<7i + <r 2 ), L = rxp. (B-4) 

The potential is a scalar under the spatial rotation. Then, V is the scalar 
functions of r, v, a\ and er 2 . 

5. Parity invariance: The strong interaction conserves parity. Thus V is invariant 
under reflection, r — > — r and v — > —v, 

V(r,v,a 1 ,a 2 ,T 1 ,T 2 ) = V(-r, -v, a 1: a 2 , n, r 2 ). (B-5) 

6. Time-reversal invariance: The strong interaction preserves time-reflection sym- 
metry under r — > r, v — > —v, ai — > —cri, which leads to 

V(r, v, cr 1; (72, Ti, r 2 ) = V(r, -v, -<ti, -<t 2 , Ti, t 2 ). (B-6) 

7. Fermi statistics: The potential is invariant under the permutation of the particle 
coordinates, 

V{r,v,cr 1 ,cr 2 ,Ti,T 2 ) = V(-r, -v,a 2 ,a 1 ,r 2 ,T 1 ) = V(r,v,<T 2 ,cr 1 ,T2,Ti), 

(B-7) 

where parity invariance was used in the second equality. 

8. Isospin invariance: The potential is invariant under the rotation in isospin space, 
which leads to two independent potentials V I=0,1 , 

V = V°(r, v, 01, 2 )P O T + V\r, v, 0i, 2 )P 1 r . (B-8) 
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9. Furthermore, V has only the terms 0""<t™ with (re, m) = (0, 0), (1, 0), (0, 1), (1, 1). 
The other higher order terms can be always reduced to the above form because 
of the property of the Pauli matrices: OiOj = 5ij + ie^kCTk- 
Then, the terms which have Pauli matrices and satisfy the above constraints are 
restricted only to the following combinations: 

c-i-o-2, (<Ti + <x 2 ) • L, (<ti • r)(cr 2 ■ r), (B-9) 
(«ri-i;)(o-2-t7), (*i-L)(*2-L). (B-10) 

It is sometimes convenient to reorganize the above 5 terms into the following hermi- 
tian operators: 

<ri-er 2 , S12 = 3(<ri • f)(<7 2 • r) - 0\ ■ a 2 , 
L S, 

P 12 = (a 1 -v)(a 2 -v), W 12 = Qi 2 -~(<t 1 -ct 2 )L 2 , (B-ll) 

where 

Qia = \ [(o-i • Z)(cr 2 • L) + (cr 2 • L)((Ti • £)] . (B-12) 

In Qi 2 the spins need to be symmetrized to make it hermitian since Lj and Lj do not 
commute with each other. Note that the term such as (S • L) 2 can be decomposed 
into Q12, S ■ L and spin-independent L 2 term. 

We decompose V° and V 1 in terms of the above operators with coefficients Vj[ 
(I = (0, 1), A = (0, a, T, LS, P, W)) which are the scalar function made of r and v 
satisfying the general constraints; 

Vi = Vi(r 2 ,v 2 ,L 2 ). (B-13) 

Note that the scalar (r • p) 2 can be written by r 2 p 2 and L 2 . 

Combining all, we arrive at the general decomposition given in 

Appendix C 

Matrix Element of the Potential 

In this appendix, we consider the partial wave decomposition of the general form 
of the ./ViV potential. At given J, there are 2 distinct states, the spin-singlet (s = 0) 
state and the spin-triplet (s = 1) state. We now consider how the five operators in 
Eq. (IB-llj) act on these states. 

The singlet state is denoted as 1 Jj, since it has s = and J = £. The fact that 
I + £ + s must be odd to satisfy fermion anti-symmetry gives 1 = for odd J and 
1 = 1 for even J. The eigenstate with J z = M can be easily obtained as 

\ 1 Jj,M) = \M,0) Jfi , (C-l) 

where we use the short-handed notation, \ J Z , s z )j tS = \ J, J z ) (g> |s, s z ). 
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The spin-triplet state is classified into 3 types: 3 7j, 3 (J ± l)j. For the first 
one, 7 = (even J) or 7 = 1 (odd J), and vice versa for other two types. By the 
Wigner-Eckart theorem, the matrix elements of the five operators do not depend on 
J z . Therefore it is enough to know eigenstates with J z = J only. Explicitly we have 



Jj,J) = -7^=j{\J-hl)j,l-^J\J,0)j, 1 } , 



\ 3 (j-i)j,j) = \j-i,i)j- 1 , 1 , 

| 3 (7 + l)j,7) 



(C-2) 
(C-3) 



:{|J-l,l)j +M 



V(7+l)(2J + 3) 
+ v / 2JTT [7(J + 1)| J + 1, -l) J+ i,i - | J, 0)7+1,!]}. (C-4) 
Using these eigenstates, it is easy to see 

(C-5) 
(C-6) 



<Ti ■ (T 2 = 2s(s + 1) - 3 = -3, 1, 1, 1 

J(J + l)-£{£ + l)-s(s + l) 



L S = 



= 0, -1, J-l, -(J + 2) 



12 



(27-l)(27 + 3) (J-l)(2J-3) (J + 2)(27 + 5) 
0, 3 , 3 , 3 ,(C7) 

for : Jj, 3 Jj, 3 (J — l)j and 3 (J+ l)j, respectively. 

For 5i2 and P12 results are more complicated due to the mixing between 3 (7— 1) j 
and 3 (7 + l)j. After a little algebra we obtain 



'12 



0, 2, 



/ 2 ( J ~ 1 ) 6^.7(7 + 1) \ 
2J + 1 ' 27+1 

6^7(7+1) 2(J + 2) 



2J + 1 
/ 2 (J-1,_, 



(C-8) 



^Pis = 0, 2p 2 , 



27 + 1 



Pj-v 



27 + 1 / 
6^/(7 + 1) 2 \ 



2J+ 1 



P+ 



where 



2 2-^i 
Pe=Pr~ l ~Pr + 

r 



6^7(7 + 1) 2 2(7 + 2) 2 
27+1 27 + 1 / 



(C-9) 



Pi 



.7 + 1 



2 2 

. J + 2 i 9 2 
1 I = \i v + , 



P- = [Pr + I 



.7 



.7 1 \ > 

Pr + * 1 ) = H V_, 



(C-10) 
(C-ll) 
(C-12) 



with £ = 7 ± 1 and ip r = 5/ (3r). 

Using these results, we obtain the potential for each channel: We have 



VfJj] = Vi{r 2 ,v 2 jJ 2 ) + U/(r 2 ,Wj, J 2 ), 7 2 = J(7 + 1) (C-13) 
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for the 1 J j state, and 

V?Jj] = Vi{r\v], J 2 ) - 3Vj(r 2 ,v 2 , j 2 ) - Vl s (r 2 ,v 2 , J 2 ) + 2y/(r 2 ,, 2 , J 2 ) 

(C-14) 



[2 J \)(2J + 6) j . 2 2 92\ Ctrl? 2 2 t2\ 2i 

~ v w( r ^J> J ) + {Vp(»" >^j> j )> u j) 



for the 3 Jj state, where 1 = 1 — 1. 

For 3 (J =F l)j, the result is more involved: 



V[ 3 (JTl)j] 



V- 



V-+ 

v ++ 



(C-15) 



where 



= V£ (r 2 , vj_ , J 2 ) - 3V< (r 2 , ,1 , J_) + (J — l)Vl s (r 2 ,vj_ , J 2 _) 
(J-l)(2J-3) j 2 2 J2 

-27TT [ 2 ^(r 2 ,-L, J 2 )+{^(r 2 ,.L,J!),^}; , 
V++ = Vq> 2 , w 2 + , J+) - 3K 7 (r 2 , v 2 + Jl) - ( J + 2)V/ S (r 2 , V 2 + , J 2 ) 
, (J + 2)(2J + 5) r, 2 , ?2 , 



(C-16) 



-^(^,«j + ,J^) 



J + 2 
~2J + 1 
3/7(7+1) 



2^(r 2 ,, 2 + ,J 2 ) + {y/(r 2 ,, 2 + ,J 2 ),, 2 + } 

7/2 2 ?2\ , ot^/„2 2 ?2\ 



(C-17) 



2(2J + 1) 

3VJ(J + 1) 
2(2J + 1) 



2Vi(r 2 ,vj + ,J^) + 2V4(r z ,vj_,J 2 _ 

+ v 2 + vUr 2 ,vj_ , J 2 ) + 4^/(r 2 , V 2 + , j 2 + )v\] , (C-18) 
2^(r 2 ,t; 2 + ,J 2 ) + 2y T / (r 2 ,.; 2 _,J 2 ) 

+^V^(r 2 ,^_, J 2 )+« 2 F/(r 2 ,t; 2 + ,J 2 )l , (C-19) 



where J± = J ± 1 and J| = J±(J± + 1). Note that V + _ = (V^) 1 with r 2 from the 
integration measure. 

Appendix D 

Heat-Kernel Representation of the Green's Function 



We define the heat kernel JC(t,x) through the initial value problem as 
d 

— fC(t, x) = V 2 IC(t, x), lim IC(t, x) = Si at {x), 
at t->o+ 



(D-l) 



where <5i at (a:) = -p- Yln&z^ e 2mn ' x l L denotes the delta function in a periodic box of 
spatial extension L. K.(t,x) is explicitly expressed as 

JC(t, x) = e tv2 5\ &t {x) = -3 Yj ex P (-i(2vr/L) 2 n 2 + 2mn ■ x/L) . (D-2) 
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For convenience, we define a modified heat kernel )C(t, x; k 2 ) as 

JC(t,x;k 2 ) = e tk2 JC(t,x), (D-3) 
which is the solution to the equation that 

TrK.(t, X] k 2 ) = (V 2 + A; 2 ) K,(t, x; k 2 ), lim K(t, x; k 2 ) = S^fx). (D-4) 
ot t->o+ 

An integration from to s leads us to 

it{s,x-k 2 ) = 5 ut (x) + f dt(V 2 + k 2 )JC(t,x;k 2 ). (D-5) 

Jo 



:k 2 ) 



This identity gives an expression for Green's function as 

G(x;k 2 ) = - (V 2 + k 2 )~ 1 IC{s,x;k 2 ) + f dtK(t,x;, 

Jo 

= - e sk2 (yi + k 2 y 1 K{s,x)+ f ' dt e tk ' X(t,x). (D-6) 

Jo 

By inserting Eq. ()D-2p . we have 

2 _ e sk2 ^ exp (-s(27r/L) 2 n 2 + 2vrm • x/L) 
{X; ' ~ ~W ^ (2Tr/L) 2 n 2 - k 2 

+ / *^j: «*(-«<■-•*>')■ (D ' 7) 

where, in the second term, we used Poisson's summation formula 

nez VP / peZ V / 

The convergence of the summations and the integration in Eq. (ID-7I) is quite good 
except at p = 0, where the integration in the second term has to be done analytically 
for t ~ 0. For this purpose, we use the formula 

1 r * «p(-£), (D-9, 



4vrr 7 (4^t) 3 / 2 ^ V 4i 
and the integration can be evaluated as 

dt e tk2 ( _ 1_ 2 

o {^tjV 2 ^ \4t X 

1 [ s dt ( tk 2 \ / a: 2 \ f 00 (it / s 2 

e — 1 exp — — — / — -r- exp 



4tt|5b| J (4vrt) 3 / 2 V / V 4t 7 7 S (4^) 3 / 2 ^ \ 4t. 

(D-10) 

Note that s dependences in Eq. (|D-6|) and Eq. (|D-7|) cancel out on the right-hand 
side, and s plays a role of the cutoff A of Eq. (D.2) in Ref. 15) as s ~ 1/A 2 . It 
controls the convergence of the summation in the first term in Eq. ()D-7p . 
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